REHO Alpraz: The effects of GABAergic agonists on regional homogeneity
Valerie Jill Sydnor
participants <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/sub_ses_passQC_finallist.csv", header = FALSE)cifti.vertex.mask <- ciftiTools::read_xifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii")glasser.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser360_regionlist.csv", header = T)schaefer.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/schaefer400_regionlist.csv", header = T)SAaxis <- read.csv("/Users/valeriesydnor/Documents/ResearchProjects/Neuron Neurodevelopmental Imaging/CorticalMap_Code/Sensorimotor_Association_Axis_AverageRanks.csv", header=T)
brains <- read.csv("/Users/valeriesydnor/Documents/ResearchProjects/Neuron Neurodevelopmental Imaging/CorticalMap_Code/brainmaps_glasser.csv", header=T)
mappys <- read.csv("/Users/valeriesydnor/Software/mappys/mappings_400.csv", header=T)Parcellate REHO Data
Parcellate vertex level REHO data into glasser360 and schaefer400 parcels for analysis
alpraz_reho_parcellate <- function(subid, sesid){
#read in xcp_abcd left and right hemisphere reho giftis
reho.lh <- readgii(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/xcpabcd/task_regress/surface/xcp_abcd/%1$s/%2$s/func/%1$s_%2$s_task-emotionid_space-fsLR_den-32k_hemi-L_desc-reho_den-91k_bold.func.gii", subid, sesid)) #32492 vertex data in $reho.lh$data$normal
reho.rh <- readgii(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/xcpabcd/task_regress/surface/xcp_abcd/%1$s/%2$s/func/%1$s_%2$s_task-emotionid_space-fsLR_den-32k_hemi-R_desc-reho_den-91k_bold.func.gii", subid, sesid)) #32492 vertex data in $reho.rh$data$normal
#combine left and right hemisphere data and save as cifti
xii.reho <- as_xifti(cortexL = reho.lh$data$normal, cortexL_mwall = cifti.vertex.mask$meta$cortex$medial_wall_mask$left, cortexR = reho.rh$data$normal, cortexR_mwall = cifti.vertex.mask$meta$cortex$medial_wall_mask$right)
dir.create(path = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/", subid, sesid), recursive = TRUE)
write_cifti(xii.reho, file.path(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii", subid, sesid)))
#parcellate vertex reho data with workbench
#GLASSER360 HCP Multimodal template parcels
command1 = sprintf("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360.pscalar.nii", subid, sesid)
ciftiTools::run_wb_cmd(command1, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
#SCHAEFER400 template parcels
command2 = sprintf("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/Schaefer2018_400Parcels_17Networks_order.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400.pscalar.nii", subid, sesid)
ciftiTools::run_wb_cmd(command2, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
}Study-Specific Cortical Mask
Create a parcellated study-specific group mask based on the fslr vertex-level group mask
#parcellate the vertex-level group mask
#GLASSER360
ciftiTools::run_wb_cmd("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_unthresholded.pscalar.nii", intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
#SCHAEFER400
ciftiTools::run_wb_cmd("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/Schaefer2018_400Parcels_17Networks_order.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_unthresholded.pscalar.nii", intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)#threshold the group mask and save out binary 0/1 mask vector
#GLASSER360
cifti.glasser.mask.unthresholded <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_unthresholded.pscalar.nii")
#information about the unthresholded mask: 191 parcels have 100% of vertices included in the group mask, 3 parcels have 50-60% of vertices included, 6 parcels have 60-70%, 5 parcels have 70-80% of vertices included, 9 parcels have 80-90% of vertices included, 18 parcelts have 90-99% of vertices included
#threshold the parcellated mask to retain only parcels where >= 50% of vertices have data for all subjects (final parcel number = 232)
glasser360.parcelmask.unthresholded <- cifti.glasser.mask.unthresholded$data
glasser360.parcelmask.thresholded <- ifelse(glasser360.parcelmask.unthresholded < 0.5, 0, glasser360.parcelmask.unthresholded)
glasser360.parcelmask.thresholded <- ifelse(glasser360.parcelmask.thresholded >= 0.5, 1, glasser360.parcelmask.thresholded)
#SCHAEFER400
cifti.schaefer.mask.unthresholded <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_unthresholded.pscalar.nii")
#information about the unthresholded mask: 178 parcels have 100% of vertices included, 5 parcels have 50-60% of vertices, 2 parcels have 60-70% of vertices, 12 parcels have 70-80%, 9 havae 80=90%
#threshold the parcellated mask to retain only parcels where >= 50% of vertices have data for all subjects (final parcel number = 227)
schaefer400.parcelmask.unthresholded <- cifti.schaefer.mask.unthresholded$data
schaefer400.parcelmask.thresholded <- ifelse(schaefer400.parcelmask.unthresholded < 0.5, 0, schaefer400.parcelmask.unthresholded)
schaefer400.parcelmask.thresholded <- ifelse(schaefer400.parcelmask.thresholded >= 0.5, 1, schaefer400.parcelmask.thresholded)#save mask as vector
write.csv(x = glasser360.parcelmask.thresholded, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_thresholded.txt", row.names = F, quote = F)
write.csv(x = schaefer400.parcelmask.thresholded, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_thresholded.txt", row.names = F, quote = F)Visualize mask
#visualize the parcellated and thresholded group mask
#GLASSER360
finalmask <- as.data.frame(glasser360.parcelmask.thresholded)
finalmask$label <- glasser.parcel.labels$label
ggseg(.data = finalmask, atlas = "glasser", mapping=aes(fill=glasser360.parcelmask.thresholded), position = c("stacked")) + theme_void() + scale_fill_gradient2(low= "black", high = "dodgerblue3", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.5) + theme(legend.position = "none")
rm(finalmask)
#SCHAEFER400
finalmask <- as.data.frame(schaefer400.parcelmask.thresholded)
finalmask$label <- schaefer.parcel.labels$label
ggseg(.data = finalmask, atlas = "schaefer17", mapping=aes(fill=schaefer400.parcelmask.thresholded), position = c("stacked")) + theme_void() + scale_fill_gradient2(low= "black", high = "mediumpurple3", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.5) + theme(legend.position = "none")
Create Masked and Parcellated REHO Data Matrices
alpraz_reho_glasser <- function(subid, sesid){
#read in parcellated reho data
cifti.reho.glasser <- read_cifti(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360.pscalar.nii", subid, sesid))
reho.glasser.unmasked <- as.array(cifti.reho.glasser$data)
reho.glasser.masked <- sweep(reho.glasser.unmasked, MARGIN = 1, glasser360.parcelmask.thresholded, `*`)
reho.glasser.masked <- as.data.frame(reho.glasser.masked)
reho.glasser.masked[reho.glasser.masked == 0] <- NA
write.table(reho.glasser.masked, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360_masked.csv", subid, sesid), row.names = F, col.names = F, quote = F)
return(reho.glasser.masked)
}reho.subxparcel.matrix.glasser <- matrix(data = NA, nrow = 84, ncol = 364)
regionheaders <- glasser.parcel.labels$orig_parcelname
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
colnames(reho.subxparcel.matrix.glasser) <- colheaders
for(row in c(1:nrow(participants))){
subid=participants[row,1]
sesid=participants[row,2]
drug=participants[row,3]
meanFD=participants[row,4]
data <- alpraz_reho_glasser(subid, sesid)
reho.subxparcel.matrix.glasser[row,] <- cbind(subid, sesid, drug, meanFD, t(data))
}write.csv(reho.subxparcel.matrix.glasser, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_glasser.csv", row.names = F, quote = F)alpraz_reho_schaefer <- function(subid, sesid){
#read in parcellated reho data
cifti.reho.schaefer <- read_cifti(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400.pscalar.nii", subid, sesid))
reho.schaefer.unmasked <- as.array(cifti.reho.schaefer$data)
reho.schaefer.masked <- sweep(reho.schaefer.unmasked, MARGIN = 1, schaefer400.parcelmask.thresholded, `*`)
reho.schaefer.masked <- as.data.frame(reho.schaefer.masked)
reho.schaefer.masked[reho.schaefer.masked == 0] <- NA
write.table(reho.schaefer.masked, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400_masked.csv", subid, sesid), row.names = F, col.names = F, quote = F)
return(reho.schaefer.masked)
}reho.subxparcel.matrix.schaefer <- matrix(data = NA, nrow = 84, ncol = 404)
regionheaders <- schaefer.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
colnames(reho.subxparcel.matrix.schaefer) <- colheaders
for(row in c(1:nrow(participants))){
subid=participants[row,1]
sesid=participants[row,2]
drug=participants[row,3]
meanFD=participants[row,4]
data <- alpraz_reho_schaefer(subid, sesid)
reho.subxparcel.matrix.schaefer[row,] <- cbind(subid, sesid, drug, meanFD, t(data))
}Visualize Cortical REHO Maps
Glasser: Average map
Visualze across-cortex regional homogeneity (Glasser)
REHO.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_glasser.csv", header=T)
REHO.glasser.map <- REHO.glasser %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
REHO.glasser.map <- colMeans(REHO.glasser.map)
REHO.glasser.map <- as.data.frame(REHO.glasser.map)
REHO.glasser.map$label <- glasser.parcel.labels$label
ggseg(.data = REHO.glasser.map, atlas = "glasser", mapping=aes(fill=REHO.glasser.map), position = c("stacked")) + theme_void() + scale_fill_gradientn(colors = viridis_pal(option="B")(10), limits = c(.45,.765))
Glasser: S-A axis correlations
Correlation with S-A Axis
##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[181:360] and SAaxis$final.rank
## S = 276670, p-value = 0.6957
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.03654347##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[1:180] and SAaxis$final.rank
## S = 278704, p-value = 0.2896
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.09959796Correlation with G1 fMRI Ratio
##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[181:360] and brains$G1.fMRI
## S = 313580, p-value = 0.05948
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1748265##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[1:180] and brains$G1.fMRI
## S = 313654, p-value = 0.01075
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.2374892Correlation with T1T2 Ratio
##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[181:360] and brains$T1T2ratio
## S = 206588, p-value = 0.01443
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2260187##
## Spearman's rank correlation rho
##
## data: REHO.glasser.map$REHO.glasser.map[1:180] and brains$T1T2ratio
## S = 188002, p-value = 0.005444
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2582577Schaefer: Average map
Visualze across-cortex regional homogeneity (Schaefer)
REHO.schaefer <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_schaefer.csv", header=T)
REHO.schaefer.map <- REHO.schaefer %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
REHO.schaefer.map <- colMeans(REHO.schaefer.map)
REHO.schaefer.map <- as.data.frame(REHO.schaefer.map)
REHO.schaefer.map$label <- schaefer.parcel.labels$label
ggseg(.data = REHO.schaefer.map, atlas = "schaefer17", mapping=aes(fill=REHO.schaefer.map), position = c("stacked")) + theme_void() + scale_fill_gradientn(colors = viridis_pal(option="B")(10), limits = c(.45,.765))
Schaefer: S-A axis correlations
Correlation with G1fMRI
##
## Spearman's rank correlation rho
##
## data: REHO.schaefer.map$REHO.schaefer.map and mappys$SA_1
## S = 2358642, p-value = 0.001502
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.2098851Correlation with T1T2 Ratio
##
## Spearman's rank correlation rho
##
## data: REHO.schaefer.map$REHO.schaefer.map and mappys$myelin
## S = 1470286, p-value = 0.0001918
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2458045Placebo versus GABA Agonist Paired T-tests
Glasser
Effect of GABA agonist on regional homogeneity per parcel (Glasser)
P-values
drug <- REHO.glasser %>% filter(drug == 0)
drug <- drug %>% select(-subid, -sesid, -drug, -meanFD)
placebo <- REHO.glasser %>% filter(drug == 1)
placebo <- placebo %>% select(-subid, -sesid, -drug, -meanFD)
#run paired t-tests on each parcel (column)
tstats <- col_t_paired(drug, placebo)num <- nrow(as.data.frame(tstats$pvalue) %>% filter(`tstats$pvalue` < 0.05))
print(sprintf("Number of significant parcels (uncorrected): %s", num))## [1] "Number of significant parcels (uncorrected): 42"P-value Histograms
#p-value histograms
ps.fdr <- p.adjust(na.omit(tstats$pvalue))
par(mfrow=c(1,2))
hist(tstats$pvalue, col = "mediumpurple4", xlab = "Paired t-test p-values", main="")
hist(ps.fdr, col = "dodgerblue4", xlab = "Paired t-test FDR-corrected p-values", main="")
T-statistic Histogram
#distribution of t-values
hist(tstats$statistic, col = "plum2", xlab = "Paired t-test t-values", main="")
T-statistic Map
REHO higher in GABA > placebo, the t-statistic is positive (yellow) REHO lower in GABA < placebo, the t-statistic is negative (purple)
tstatistics <- tstats$statistic
tstatistics <- as.data.frame(tstatistics)
tstatistics$label <- glasser.parcel.labels$label
ggseg(.data = tstatistics, atlas = "glasser", mapping=aes(fill=tstatistics), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.88)
Correlation with S-A Axis
Correlation with S-A Axis (Glasser)
## Warning in cor.test.default(tstatistics$tstatistics[181:360],
## SAaxis$final.rank, : Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: tstatistics$tstatistics[181:360] and SAaxis$final.rank
## S = 373871, p-value = 7.575e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4007081## Warning in cor.test.default(tstatistics$tstatistics[1:180], SAaxis$final.rank, :
## Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: tstatistics$tstatistics[1:180] and SAaxis$final.rank
## S = 378506, p-value = 2.119e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4933579Correlation with T1T2 Ratio (Glasser)
##
## Spearman's rank correlation rho
##
## data: tstatistics$tstatistics[181:360] and brains$T1T2ratio
## S = 199462, p-value = 0.006101
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2527162##
## Spearman's rank correlation rho
##
## data: tstatistics$tstatistics[1:180] and brains$T1T2ratio
## S = 171112, p-value = 0.0004238
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3248954Schaefer
Effect of GABA agonist on regional homogeneity per parcel (Schaefer)
P-values
drug <- REHO.schaefer %>% filter(drug == 0)
drug <- drug %>% select(-subid, -sesid, -drug, -meanFD)
placebo <- REHO.schaefer %>% filter(drug == 1)
placebo <- placebo %>% select(-subid, -sesid, -drug, -meanFD)
#run paired t-tests on each parcel (column)
tstats <- col_t_paired(drug, placebo)## Warning: col_t_paired: 173 of the columns had less than 2 paired observations.
## First occurrence at column 23num <- nrow(as.data.frame(tstats$pvalue) %>% filter(`tstats$pvalue` < 0.05))
print(sprintf("Number of significant parcels (uncorrected): %s", num))## [1] "Number of significant parcels (uncorrected): 54"P-value Histograms
#p-value histograms
ps.fdr <- p.adjust(na.omit(tstats$pvalue))
par(mfrow=c(1,2))
hist(tstats$pvalue, col = "mediumpurple4", xlab = "Paired t-test p-values", main="")
hist(ps.fdr, col = "dodgerblue4", xlab = "Paired t-test FDR-corrected p-values", main="")
T-statistic Histogram
#distribution of t-values
hist(tstats$statistic, col = "plum2", xlab = "Paired t-test t-values", main="")
T-statistic Map
REHO higher in GABA > placebo, the t-statistic is positive (yellow) REHO lower in GABA < placebo, the t-statistic is negative (purple)
tstatistics <- tstats$statistic
tstatistics <- as.data.frame(tstatistics)
tstatistics$label <- schaefer.parcel.labels$label
ggseg(.data = tstatistics, atlas = "schaefer17", mapping=aes(fill=tstatistics), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.88)
Correlation with S-A axis
Correlation with G1 fMRI (Schaefer)
##
## Spearman's rank correlation rho
##
## data: tstats$statistic and mappys$SA_1
## S = 2921866, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4987956Correlation with T1T2 Ratio (Schaefer)
##
## Spearman's rank correlation rho
##
## data: tstats$statistic and mappys$myelin
## S = 1134078, p-value = 7.194e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4182652Correlation with PVALB (Schaefer)
##
## Spearman's rank correlation rho
##
## data: tstats$statistic and mappys$PVALB
## S = 1305640, p-value = 4.154e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3302611Correlation with SST (Schaefer)
##
## Spearman's rank correlation rho
##
## data: tstats$statistic and mappys$SST
## S = 2558606, p-value = 1.776e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.3124583REHO PCA and Logistic Regressions
Glasser
REHO PCA and PC scores logistic regression for drug/placebo prediction (Glasser)
PCA Variance Explained
inputmatrix <- REHO.glasser
inputmatrix <- inputmatrix %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
inputmatrix <- inputmatrix %>% mutate_if(is.character,as.numeric)
inputmatrix <- as.data.frame(inputmatrix) #convert matrix to df
inputmatrix[inputmatrix == 0] <- NA
inputmatrix.cleaned <- inputmatrix[, unlist(lapply(inputmatrix, function(x) !all(is.na(x))))] #remove all columns (parcels) with null data (all NAs) #84x232
reho.pca.glasser <- prcomp(inputmatrix.cleaned, scale. = TRUE, center = TRUE)
print(summary(reho.pca.glasser)$importance[2,1:10]*100)## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 25.324 5.610 5.445 4.826 3.371 3.138 2.807 2.740 2.340 2.259pve <- 100*(reho.pca.glasser$sdev)^2/sum ((reho.pca.glasser$sdev)^2)
plot(pve, pch=16,
xlab="Principal Components",
ylab="Prop. of variance explained")
Logistic Regression
reho.pca.glasser.scores <- as.data.frame(reho.pca.glasser$x[,1:10])
reho.pca.glasser.scores$subid <- participants$V1
reho.pca.glasser.scores$sesid <- participants$V2
reho.pca.glasser.scores$drug <- participants$V3
reho.pca.glasser.scores$meanFD <- participants$V4 #match subids to scores
summary(glm(as.factor(reho.pca.glasser.scores$drug) ~ reho.pca.glasser.scores$PC1 + reho.pca.glasser.scores$meanFD, family=binomial(logit)))##
## Call:
## glm(formula = as.factor(reho.pca.glasser.scores$drug) ~ reho.pca.glasser.scores$PC1 +
## reho.pca.glasser.scores$meanFD, family = binomial(logit))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.42708 -1.15218 0.06032 1.15661 1.40665
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.65963 0.60666 1.087 0.277
## reho.pca.glasser.scores$PC1 0.02048 0.02952 0.694 0.488
## reho.pca.glasser.scores$meanFD -3.30648 2.84188 -1.163 0.245
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.45 on 83 degrees of freedom
## Residual deviance: 114.26 on 81 degrees of freedom
## AIC: 120.26
##
## Number of Fisher Scoring iterations: 4Schaefer
REHO PCA and PC scores logistic regression for drug/placebo prediction (Schaefer)
PCA Variance Explained
inputmatrix <- REHO.schaefer
inputmatrix <- inputmatrix %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
inputmatrix <- inputmatrix %>% mutate_if(is.character,as.numeric)
inputmatrix <- as.data.frame(inputmatrix) #convert matrix to df
inputmatrix[inputmatrix == 0] <- NA
inputmatrix.cleaned <- inputmatrix[, unlist(lapply(inputmatrix, function(x) !all(is.na(x))))] #remove all columns (parcels) with null data (all NAs) #84x232
reho.pca.schaefer <- prcomp(inputmatrix.cleaned, scale. = TRUE, center = TRUE)
print(summary(reho.pca.schaefer)$importance[2,1:10]*100)## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 26.054 6.463 5.366 4.711 3.922 3.170 2.893 2.654 2.543 2.226pve <- 100*(reho.pca.schaefer$sdev)^2/sum ((reho.pca.schaefer$sdev)^2)
plot(pve, pch=16,
xlab="Principal Components",
ylab="Prop. of variance explained")
Logistic Regression
reho.pca.schaefer.scores <- as.data.frame(reho.pca.schaefer$x[,1:10])
reho.pca.schaefer.scores$subid <- participants$V1
reho.pca.schaefer.scores$sesid <- participants$V2
reho.pca.schaefer.scores$drug <- participants$V3
reho.pca.schaefer.scores$meanFD <- participants$V4 #match subids to scores
summary(glm(as.factor(reho.pca.schaefer.scores$drug) ~ reho.pca.schaefer.scores$PC1 + reho.pca.schaefer.scores$meanFD, family=binomial(logit)))##
## Call:
## glm(formula = as.factor(reho.pca.schaefer.scores$drug) ~ reho.pca.schaefer.scores$PC1 +
## reho.pca.schaefer.scores$meanFD, family = binomial(logit))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.46086 -1.14744 0.04421 1.14899 1.38829
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.63966 0.60842 1.051 0.293
## reho.pca.schaefer.scores$PC1 0.02505 0.02963 0.845 0.398
## reho.pca.schaefer.scores$meanFD -3.20496 2.84927 -1.125 0.261
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.45 on 83 degrees of freedom
## Residual deviance: 114.02 on 81 degrees of freedom
## AIC: 120.02
##
## Number of Fisher Scoring iterations: 4SVM
Functions
SVM_2class <- function(df,folds,feature_selection = F,feature_proportion = .1,num_repetitions = 100){
#SVM 2-class classifier
# folds: number of folds for cv.
# feature_selection: optional flag for data-driven feature selection. Selects the proportion of features indicated by feature_proportion. Not used by default.
# num_repetitions: How many times to you want to repeat the cv process using different random splits of the data. This is just extra caution against a randomly good or bad split.
cat('\nRunning SVM models.....')
# Set up folds for CV
if (folds == "LOO") {
# This is for leave-one-out cross-validation
num_folds = length(unique(df$subid))
num_repetitions <- 1
} else {
num_folds = folds
num_repetitions <- num_repetitions
}
svm_output <- vector(mode = "list",length = num_repetitions) #set up output object.
# Create the folds
unique_IDs <- unique(df$subid)
subid_folds <- replicate(num_repetitions,sample(unique_IDs,size = length(unique_IDs))) # create sets of random draws.
foldIdxs <- data.frame(subid=unique_IDs)
foldIdxs$foldID <- row_number(foldIdxs$subid)
foldIdxs$foldID <- ntile(foldIdxs$foldID,num_folds)
# cat('Sending data to CV')
for (r in 1:num_repetitions) {
foldIdxs$subid <- subid_folds[,r] # Grab a random split for folds.
fold_output<-vector(mode = "list",length = num_folds)
cat(sprintf('\nrepetition %d.... ',r))
for (fold in 1:num_folds) {
# cat(sprintf('\nfold %d.... ',fold))
trainingIDs <- as.matrix(foldIdxs %>% filter(foldID != fold) %>% select(subid))
trainingIndex <- df$subid %in% trainingIDs # indices for training subs
trainingData <- df[trainingIndex, 3:dim(df)[2] ] %>% arrange(drug) #this is important because libsvm automatically makes the first observation class 1, so drug must be first every time. Placebo will be class -1.
testData <- df[!trainingIndex, 4:dim(df)[2]] # test data. Take columns 4:end (Removes subid sesid drug).
testLabels <- data.frame(df[!trainingIndex,c(1:3) ])
# svm
x <- as.matrix(trainingData[, 2:dim(trainingData)[2]])
y <- as.factor(as.matrix(trainingData[,1]))
svm.model <- svm(x =x, y = y, cost = 1, kernel = "linear",type = "C-classification",scale = F)
svm.pred <- predict(svm.model, as.matrix(testData))
w <- t(svm.model$coefs) %*% svm.model$SV #calculate feature weights.
# num_features <- dim(x)[2]
decisionValues <- w %*% t(testData)-svm.model$rho # Get decision valus.
distance <- decisionValues/norm(w) #calculate distance from the classification hyperplane
# just adding the results to the dataframe.
testLabels$decisionValues <- t(decisionValues)
testLabels$distance <- t(distance)
testLabels$model.pred = svm.pred
fold_output[[fold]]<-testLabels
}
svm_output[[r]] <- data.table::rbindlist(fold_output) # saving output for the repetition.
cat('complete\n')
}
# output_list <- apply(subid_folds,CV_function,df=df,num_folds = num_folds,MARGIN = 2)
# Now train a model using all the data. This is for the estimation of the feature weights using the most possible data.
final_data<-df[, 3:dim(df)[2] ]
x <- as.matrix(final_data[, 2:dim(final_data)[2]])
y <- as.factor(as.matrix(final_data[,1]))
svm.model <- svm(x = x, y = y,
cost = 1, kernel = "linear",type = "C-classification",scale = F)
w <- t(svm.model$coefs) %*% svm.model$SV
svm_results <- list(svm_output,w,svm.model)
cat('\nFinished SVM models\n')
return(svm_results)
}AUC <- function(DecisionValues, labels){
# This function creates an ROC curve and then calculates AUC.
# Decision values is nx1
# Labels is nx1
# N.B.
# Drug (class 0) is assigned as +1 by libsvm and placebo (class 1) is assigned as 0 by libsvm default.
# Adjust the true labels and predicted labels to match that here so that the decision values make sense.
labels <- labels*-1+1
P <- sum(labels == 1)
N <- sum(labels == 0)
Sorted_DecisionValues <- sort(unique(DecisionValues), decreasing = FALSE)
numDecisionValues <- length(Sorted_DecisionValues)
TP_Array <- vector(mode = "numeric",length = numDecisionValues)
FP_Array <- vector(mode = "numeric",length = numDecisionValues)
Accuracy_Array = vector(mode = "numeric",length = numDecisionValues)
for (i in 1:numDecisionValues){
thisCutoff <- Sorted_DecisionValues[i]
thisPredictedLabels <- as.numeric(DecisionValues>thisCutoff)
detections <- thisPredictedLabels==1
TP <- sum(labels[detections] == thisPredictedLabels[detections])
TPR <- TP/P
FP <- sum(labels[detections]!=thisPredictedLabels[detections])
FPR <- FP/N
TP_Array[i] <- TPR
FP_Array[i] <- FPR
Accuracy_Array[i] = (TP + N - FP) / (P + N)
}
ROC_output <- data.frame(TPR =TP_Array,FPR=FP_Array,Accuracy = Accuracy_Array)
ROC_output <- ROC_output%>%arrange(TPR,FPR)
#AUC
dFPR <- c(0,diff(ROC_output$FPR))
dTPR <- c(0,diff(ROC_output$TPR))
AUC <- sum(ROC_output$TPR * dFPR) + sum(dTPR * dFPR)/2
return(AUC)
}Glasser
Drug v. Placebo SVM (Glasser)
regionheaders <- glasser.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
df <- REHO.glasser
colnames(df) <- colheaders
df <- df %>% dplyr::select(-meanFD)
df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]model_results <- SVM_2class(df, "LOO") # model_results contains [[1]] svm_output table of subid, sesid, drug, decisionValues, distance, and model.pred, [[2]] weights for each region from the model run with all sessions and [[3]] the model itself
saveRDS(model_results, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.glasser.rds")model_results <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.glasser.rds")
prediction_output <- model_results[[1]]W <- model_results[[2]] #weights for each parcel
svm.model <- model_results[[3]] #the model run
num_features <- sum(!is.na(W[1,]))Model accuracy
accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy.
accuracies <- sapply(prediction_output,accuracy_fun) #Apply the function to each repetition of the cross-validation. For LOO there is only 1 repetition, and accuracies = accuracy
accuracy <- mean(accuracies)
accuracy## [1] 0.6190476Number of correct classifications
num_obs <- length(prediction_output[[1]]$model.pred)
num_correct <- round(accuracy*num_obs)
num_correct## [1] 52Binomial p-value significance
b <- binom.test(num_correct,num_obs,.5) #This is a binomial test. The p-value is not used (permutation test is used instead), but this function summarizes the data nicely.
b##
## Exact binomial test
##
## data: num_correct and num_obs
## number of successes = 52, number of trials = 84, p-value = 0.03753
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.5065562 0.7228868
## sample estimates:
## probability of success
## 0.6190476Permutation p-value significance
num_permutations = 1000
nw=1
folds="LOO"
perm_acc <- matrix(nrow = num_permutations)
perm_W <- matrix(nrow = num_permutations,ncol = num_features)
perm_list <- list()
perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
.combine=c,
.export = c("featureExtraction","SVM_2class"),
.packages = c("dplyr","e1071","matrixTests")) %do% { # This must be `dopar` to be parallel.
perm_result_list=list()
for (p in 1:perm_chunk){
# thisPerm <- perms[p,]
cat(sprintf("\nPermutation %d",p))
thisDF <- df
# thisDF$drug <- df$drug[thisPerm]
permuted <- df %>% select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
thisDF$drug <- permuted$perm_drug
perm_pred_result <- SVM_2class(thisDF,folds,
feature_selection = feature_selection,
feature_proportion = feature_proportion,
num_repetitions = 1)
perm_result_list[[p]] = list(perm_pred_result)
}
perm_result_list
}
#organizing the output
pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])# Calculating permutation accuracies and AUCs
perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))
perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)
perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))
W_test <- sapply(perm_Ws, FUN = function(x) {
abs(x)>abs(W)
})
W_sig <- rowMeans(W_test)
b$perm_p <- perm_p
b$perm_W_sig <- W_sig
b$perm_accs = perm_acc_distribution
b$perm_aucs = perm_auc_distribution
saveRDS(b, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.glasser.rds")
saveRDS(W, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.glasser.rds")b <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.glasser.rds")
W <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.glasser.rds")
print("How many times is the permuted accuracy greater than the actual accuracy?")## [1] "How many times is the permuted accuracy greater than the actual accuracy?"## [1] 1## [1] "Overall Accuracy: 0.619; p = 0.00100\n\n"P-value distribution
Model weights map
SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$label <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"
ggseg(.data = SVM.weights.glasser, atlas = "glasser", mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)
Model weights S-A axis correlation
SAaxis$label <- glasser.parcel.labels$label[181:360]
SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)##
## Spearman's rank correlation rho
##
## data: SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 343225, p-value = 0.00178
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.2858925Model accuracy compared to Catch22 Features
catch22_SVM_accuracy <- function(feature){
catchnames <- (c("DN_HistogramMode_5","DN_HistogramMode_10","CO_f1ecac","CO_FirstMin_ac",
"CO_HistogramAMI_even_2_5", "CO_trev_1_num", "MD_hrv_classic_pnn40",
"SB_BinaryStats_mean_longstretch1", "SB_TransitionMatrix_3ac_sumdiagcov",
"PD_PeriodicityWang_th0_01", "CO_Embed2_Dist_tau_d_expfit_meandiff",
"IN_AutoMutualInfoStats_40_gaussian_fmmi", "FC_LocalSimple_mean1_tauresrat",
"DN_OutlierInclude_p_001_mdrmd", "DN_OutlierInclude_n_001_mdrmd",
"SP_Summaries_welch_rect_area_5_1", "SB_BinaryStats_diff_longstretch0",
"SB_MotifThree_quantile_hh", "SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1",
"SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1", "SP_Summaries_welch_rect_centroid",
"FC_LocalSimple_mean3_stderr"))
#prepare data for SVM
catchname <- catchnames[feature]
catchdata <- read.csv(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/catch%s_sesxparcel_matrix.csv", feature))
df1 <- catchdata[,1:3]
df2 <- catchdata[,5:364]
df2 <- df2 %>% mutate_if(is.character,as.numeric)
df <- cbind(df1, df2)
df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]
#Run SVM
catchmodel_results <- SVM_2class(df, "LOO") # model_results contains [[1]] svm_output table of subid, sesid, drug, decisionValues, distance, and model.pred, [[2]] weights for each region from the model run with all sessions and [[3]] the model itself
saveRDS(catchmodel_results, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch%s.Rdata", feature))
catchprediction_output <- catchmodel_results[[1]]
catchprediction_weights <- catchmodel_results[[2]] #weights for each parcel
#Calculate accuracy and binomial p-value
accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy.
accuracy <- sapply(catchprediction_output,accuracy_fun)
num_obs <- length(catchprediction_output[[1]]$model.pred)
num_correct <- round(accuracy*num_obs)
b <- binom.test(num_correct,num_obs,.5) #This is a binomial test
binomial_pvalue <- b$p.value
#Calculate permutation-based p-values
num_permutations = 1000
nw=1
folds="LOO"
num_features = 232
perm_acc <- matrix(nrow = num_permutations)
perm_W <- matrix(nrow = num_permutations,ncol = num_features)
perm_list <- list()
perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
.combine=c,
.export = c("featureExtraction","SVM_2class"),
.packages = c("dplyr","e1071","matrixTests")) %do% { # This must be `dopar` to be parallel.
perm_result_list=list()
for (p in 1:perm_chunk){
# thisPerm <- perms[p,]
cat(sprintf("\nPermutation %d",p))
thisDF <- df
# thisDF$drug <- df$drug[thisPerm]
permuted <- df %>% select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
thisDF$drug <- permuted$perm_drug
perm_pred_result <- SVM_2class(thisDF,folds,
feature_selection = feature_selection,
feature_proportion = feature_proportion,
num_repetitions = 1)
perm_result_list[[p]] = list(perm_pred_result)
}
perm_result_list
}
pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])
perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))
saveRDS(perm_acc_distribution, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_permutedaccuracies_catch%s.Rdata", feature))
perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)
perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))
compiled.output <- cbind(catchname, accuracy, binomial_pvalue, perm_p)
return(compiled.output)
}catch22.SVM.accuracy <- matrix(NA, ncol = 4, nrow = 22)
for(feature in c(1:22)){
thisfeature.accuracy <- catch22_SVM_accuracy(feature)
catch22.SVM.accuracy[feature,] <- thisfeature.accuracy
}
colnames(catch22.SVM.accuracy) <- c("Catch Feature", "SVM Accuracy", "Binomial pval", "Permutation pval")write.csv(catch22.SVM.accuracy, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/Catch22_SVM_Accuracies.csv", quote = F, row.names = F)
catch22.SVM.accuracycatch22.SVM.accuracy <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/Catch22_SVM_Accuracies.csv")hist(as.numeric(catch22.SVM.accuracy[,2]), col="grey", breaks=12, xlim = c(0.4,0.65))
abline(v = accuracy, col="red", lwd=3, lty=2)
Best Catch22 Features
Catch9
catch9.modelresults <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch9.Rdata")
prediction_output <- catch9.modelresults[[1]]
W <- catch9.modelresults[[2]]SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$origlabel <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label <- NA
for(weight in c(1:232)){
newlabel <- substr(SVM.weights.glasser$origlabel[weight],1,nchar(SVM.weights.glasser$origlabel[weight])-4)
if(substring(newlabel, 1, 1) == "R"){
newlabelhemi <- gsub("R_", "rh_R_", newlabel)}
else
newlabelhemi <- gsub("L_", "lh_L_", newlabel)
SVM.weights.glasser[weight,3] <- newlabelhemi
}SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)##
## Spearman's rank correlation rho
##
## data: SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 269943, p-value = 0.1958
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1225927SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9.46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"
ggseg(.data = SVM.weights.glasser, atlas = "glasser", mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)## merging atlas and data by 'label'## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label V1 origlabel
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_9.46d -0.0154 L_9.46d_ROI
## 2 <NA> <NA> <NA> <NA> <NA> lh_L_10pp -0.0431 L_10pp_ROI
Catch13
catch13.modelresults <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch13.Rdata")
prediction_output <- catch13.modelresults[[1]]
W <- catch13.modelresults[[2]]SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$origlabel <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label <- NA
for(weight in c(1:232)){
newlabel <- substr(SVM.weights.glasser$origlabel[weight],1,nchar(SVM.weights.glasser$origlabel[weight])-4)
if(substring(newlabel, 1, 1) == "R"){
newlabelhemi <- gsub("R_", "rh_R_", newlabel)}
else
newlabelhemi <- gsub("L_", "lh_L_", newlabel)
SVM.weights.glasser[weight,3] <- newlabelhemi
}SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)##
## Spearman's rank correlation rho
##
## data: SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 214410, p-value = 0.2533
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1083493SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"
ggseg(.data = SVM.weights.glasser, atlas = "glasser", mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)## merging atlas and data by 'label'## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label V1 origlabel
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_10pp -0.499 L_10pp_ROI
Schaefer
Drug v. Placebo SVM (Schaefer)
regionheaders <- schaefer.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
df <- REHO.schaefer
colnames(df) <- colheaders
df <- df %>% select(-meanFD)
df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]model_results <- SVM_2class(df, "LOO")
saveRDS(model_results, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.schaefer.rds")model_results <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.schaefer.rds")
prediction_output <- model_results[[1]]Model accuracy
accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy.
accuracies <- sapply(prediction_output,accuracy_fun) #Apply the function to each repetition of the cross-validation.
accuracy <- mean(accuracies)
accuracy## [1] 0.6309524Number of correct classifications
num_obs <- length(prediction_output[[1]]$model.pred)
num_correct <- round(accuracy*num_obs)
num_correct## [1] 53Binomial p-value significance
b <- binom.test(num_correct,num_obs,.5) #This is a binomial test. The p-value is not used (permutation test is used instead), but this function summarizes the data nicely.
b##
## Exact binomial test
##
## data: num_correct and num_obs
## number of successes = 53, number of trials = 84, p-value = 0.02138
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.5186641 0.7337192
## sample estimates:
## probability of success
## 0.6309524Permutation p-value significance
num_permutations = 1000
nw=1
folds="LOO"
perm_acc <- matrix(nrow = num_permutations)
perm_W <- matrix(nrow = num_permutations,ncol = num_features)
perm_list <- list()
perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
.combine=c,
.export = c("featureExtraction","SVM_2class"),
.packages = c("dplyr","e1071","matrixTests")) %do% { # This must be `dopar` to be parallel.
perm_result_list=list()
for (p in 1:perm_chunk){
# thisPerm <- perms[p,]
cat(sprintf("\nPermutation %d",p))
thisDF <- df
# thisDF$drug <- df$drug[thisPerm]
permuted <- df %>% dplyr::select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
thisDF$drug <- permuted$perm_drug
perm_pred_result <- SVM_2class(thisDF,folds,
feature_selection = feature_selection,
feature_proportion = feature_proportion,
num_repetitions = 1)
perm_result_list[[p]] = list(perm_pred_result)
}
perm_result_list
}
#organizing the output
pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])# Calculating permutation accuracies and AUCs
perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))
perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)
perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))
W_test <- sapply(perm_Ws, FUN = function(x) {
abs(x)>abs(W)
})
W_sig <- rowMeans(W_test)
b$perm_p <- perm_p
b$perm_W_sig <- W_sig
b$perm_accs = perm_acc_distribution
b$perm_aucs = perm_auc_distribution
saveRDS(b, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.schaefer.rds")
saveRDS(W, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.schaefer.rds")b <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.schaefer.rds")
W <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.schaefer.rds")
print("How many times is the permuted accuracy greater than the actual accuracy?")## [1] "How many times is the permuted accuracy greater than the actual accuracy?"## [1] 2## [1] "Overall Accuracy: 0.631; p = 0.00200\n\n"P-value distribution
Model weights S-A axis correlation
SVM.weights.schaefer <- as.data.frame(t(W))
SVM.weights.schaefer$label <- row.names(SVM.weights.schaefer)
mappys$label <- schaefer.parcel.labels$label
SVM.weights.schaefer.SA <- merge(SVM.weights.schaefer, mappys, by="label")
cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$SA_1, method=c("spearman"), exact=F)##
## Spearman's rank correlation rho
##
## data: SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$SA_1
## S = 2737344, p-value = 2.493e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4041435Model weights inhibitory interneuron correlation
##
## Spearman's rank correlation rho
##
## data: SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$PVALB
## S = 1504824, p-value = 0.0005336
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.228088##
## Spearman's rank correlation rho
##
## data: SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$SST
## S = 2286410, p-value = 0.009073
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1728331##
## Spearman's rank correlation rho
##
## data: SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$GABRA3
## S = 2280092, p-value = 0.01048
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1695922Model weights NMDA receptor subtype distribution
##
## Spearman's rank correlation rho
##
## data: SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$GRIN2B
## S = 2187336, p-value = 0.0665
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1220123